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A FINITE ELEMENT - BOUNDARY INTEGRAL METHOD FOR 


SCATTERING AND RADIATION BY TWO- AND 
THREE-DIMENSIONAL STRUCTURES 

Jian-Ming Jin, John L. Volakis , and Jeffery D. Collins 
Radiation Laboratory 

Department of Electrical Engineering and Computer Science 
The University of Michigan 
Ann Arbor , Michigan 48109-2122 

ABSTRACT 

This paper presents a review of a hybrid finite element - boundary integral formu- 
lation for scattering and radiation by two- and three-dimensional composite structures. 
In contrast to other hybrid techniques involving the finite element method, the proposed 
one is in principle exact and can be implemented using a low 0(N) storage. This is of 
particular importance for large scale applications and is a characteristic of the boundary 
chosen to terminate the finite element mesh, usually as close to the structure as possible. 
A certain class of these boundaries lead to convolutional boundary integrals which can 
be evaluated via the fast Fourier transform (FFT) without a need to generate a matrix, 
thus, retaining the O(N) storage requirement. The paper begins with a general descrip- 
tion of the method. A number of two- and three-dimensional applications are then given, 
including numerical computations which demonstrate the method’s accuracy, efficiency 
and capability. 



I. INTRODUCTION 


Integral equation methods and to a lesser degree differential equation methods are 
commonly used for electromagnetic scattering and radiation computations. Among 
the various integral equation approaches, surface formulations (see, for example, [1]- 
[4]) are more attractive for perfectly conducting, impedance or layered material struc- 
tures, whereas volume formulations (see, for example, [5]- [8]) are particularly suited 
for modeling inhomogeneous scatterers. Both have been applied to a variety of two- 
and three-dimensional problems and are often referred to as exact techniques because 
they guarantee convergence for sufficiently dense discretizations. However, they have 
the disadvantage of being difficult to implement for complex objects and also result in 
full matrices whose treatment requires a large memory. This is particularly true for 
three-dimensional applications and because of it, differential equation approaches are 
becoming more popular. 

Differential equation methods can be subdivided into finite element and finite differ- 
ence methods. In contrast to the integral equation approaches, they lead to relatively 
simple formulations and are thus attractive for simulating complex penetrable structures. 
More importantly, they are associated with sparse, banded matrices which can be effi- 
ciently solved and stored. They do not, however, incorporate the Sommerfeld radiation 
condition and this requires that the domain of discretization be extended far from the 
scatterer where the radiation condition can be imposed [9]. This is a major disadvantage 
of the differential equation methods and recent efforts have concentrated on the use of 
absorbing boundary conditions to reduce the discretization region outside the scatterer 
(see, for example, [10], [11]). Unfortunately, the accuracy of the absorbing boundary 
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conditions is dependent upon the composition and shape of the scatterer, leading to 
results of unpredictable accuracy. 

To eliminate the disadvantages of the integral and differentia] equation methods while 
retaining their advantages, various hybrid methodologies have been developed (see, for 
example, [12]- [17]), The general principle of hybrid techniques is to introduce a fictitious 
boundary enclosing the scatterer. Interior to the boundary, the finite element or finite 
difference method is used to formulate the fields whereas in the exterior region the 
fields are represented by an eigenfunction expansion [12]-[14] or a boundary integral [15]- 
[17]. The last is an exact representation, but in practice the eigenfunction expansion is 
approximate since the infinite series must be truncated. 

In this paper, we describe a hybrid technique which combines the finite element 
and boundary integral methods. The paper is essentially a review of our recent work 
pertaining to the development of a finite element - boundary integral method [18]-[24], 
In the next section we present the general formulation without reference to any specific 
geometry or application. A number of two- and three-dimensional applications are then 
considered to demonstrate the accuracy, efficiency and capability of the method. Of 
particular concern in these applications is the choice of the fictitious boundary enclosing 
the structure so that the resulting boundary integrals are convolutions. They can then be 
evaluated via the fast Fourier transform (FFT) in conjunction with an iterative solution 
approach such as the conjugate gradient (CG) or biconjugate gradient (BiCG) method. 
In this manner, the generation of a partly full matrix is avoided and the solution requires 
only O(N) storage without compromise in accuracy. Also, provided the enclosure fits 
tightly over the structure, the elements comprising the mesh outside the domain of the 
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structure are kept to a minimum. 


II. GENERAL FORMULATION 

Consider the scattering/radiation problem illustrated in Figure 1. We are interested 
in the computation of the fields generated by the external sources ( J ext , M cx *) (for the 
scattering case) or by the internal sources (J m< , M tn *) (for the radiation case) in the 
presence of a three-dimensional structure immersed in an infinite, homogeneous medium. 
In the following we describe a numerical procedure for determining the field everywhere 
by combining the finite element and boundary integral methods. 

1. Decoupling and Coupling of Exterior and Interior Fields 

To combine the finite element and boundary integral methods, it is necessary that the 
three-dimensional structure be enclosed in a fictitious surface denoted by S. Within S 
the finite element method is employed to formulate the fields, whereas outside S the fields 
are represented by the radiation of the extenal sources and a set of equivalent electric 
and magnetic currents placed on S. This permits the decoupling of the fields interior 
and exterior to S which can later be coupled by enforcing tangential field continuity on S 
leading to a system of equations for the solution of the equivalent electric and magnetic 
currents. 

By making use of the free space Green’s function, the fields in the region exterior to 
S (this region will be hereon denoted as V^) are represented as 

E (or H) = E inc (or H ,nc ) + (E + x n, H + x n) (1) 

where (E ,nc ,H‘ nc ) denote the incident fields radiated by the impressed external sources 
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(J ex *, M ext ) and denote the boundary/surface operators (integral or integro-differential) 
yielding the field radiated by the equivalent electric (H + x n) and magnetic (E + x h) 
currents. The subscript e corresponds to the operator associated with the E field formu- 
lation and likewise the subscript h is associated with the H field representation. Also, 
(E + ,H+) are the fields on S as one approaches from and n denotes the unit vector 
normal to S and pointing toward Let us now consider the region interior to 5, 
denoted as V, When this region is inhomogeneous, its associated Green’s function is 
usually not available and we cannot, therefore, formulate the fields in terms of bound- 
ary integrals as was done for However, by using the differential form of Maxwell’s 
equations, we can find a relation between the interior and the tangential boundary fields. 
This can be of the form 

h t h (E or H, E" x n, H“ xn)=0 (2) 

where f Ct k denote the appropriate operators and (E“,H~) are the fields on S as one 
approaches from V. We note that, since / e ^ do not involve any Green’s function, their 
discretization via the finite element method results in sparse and banded matrices. Con- 
sequently, the implementation of (2) is associated with 0(N) memory demand which is 
one of the primary advantages of the finite element method. 

To solve for (E± x n, H ± x n) we must couple the fields given in (1) and (2) by 
enforcing continuity of the tangential electric and magnetic fields on 5. The continuity 
condition demands that 

E + x n = E~ x n , H + x n = H“ X n (3) 

which together with (1) and (2) imply a system for the solution of the interior and 
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boundary fields. 


2. Formulation of Functionals and Weighted Residual Equations 

To formulate the fields in V, we begin with the vector wave equation 


V x x E j - jfcge r E = -jk 0 Z 0 3 int + V x (4) 

where ko = 2tt/\ is the free space wavenumber and Zo is the free space intrinsic 
impedance. A traditional approach for solving this is to consider the functional 


-F(E) = | JJJ v i(VxE).(VxE)-^ r E.E dV 


+ JJJ r E • [jk 0 Z 0 J int - V x 

+jk 0 Z 0 jj)^ E • (H x n ) dS 


dV 


(5) 


which can be easily shown to be stationary with respect to the solution of (4) with 
(H x n) being considered as one of the sources for E. Thus E can be found by enforcing 


SF( E) = 0 


( 6 ) 


where 6F( E) denotes the first order variation of F about E. 

An alternative approach for solving (4) is to employ the method of weighted residuals. 
This is an approach often presented in graduate electromagnetics texts and it would, 
therefore, be instructive to employ it here as well. Based on the weighted residual 
method we demand that 


<R,T) = JJJ v Vx^VxEj.T-4E.T dV 




jk 0 z 0 r nt - v x 




dV = 0 


(7) 
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where R denotes the residual of (4) and T is a testing or weighting function chosen to 
satisfy the required boundary conditions associated with E. Recalling the identity 


V/ir 

= JJJ v -j-(V x E) . (V x T) dV + jk oZ 0 jj s T • (H x n) dS (8) 


V x —V x E ) • T dV 

V V^r 


(7) can be written as 


<R,T> = 


JJJ v [^-(V x E) • (V x T) - k 2 0 e r E • t] dV 
+;Jt 0 Z 0 j/jf T • (H x n) dS = 0 


( 9 ) 


which is equivalent to (5) and (6) provided T is chosen to be the same as the expansion 
basis, implying an application of Galerkin’s technique. 

In some cases, it may be advantageous to work with magnetic fields rather than 
electric fields. It will then be necessary to use the dual of (5) and (9) given by 


F( H) = ^ JJJ v ^-(V x H) • (V x H) - k 2 fi T Il • H dV 


H* 


jk 0 Y 0 M int + V 


x 


dV 


and 


-III 

—jkoYo jj)$ H • (E x n) dS 


(R, T) = JJJ v [i(V x H) . (V x T) - kforK * t] dV 

- JIIv T * ^° y ° M,nt + v x (f J ‘ n< )] 


( 10 ) 



respectively, where Yq = 1/Zo- 


3. Finite Element Discretization 

To discretize V , we subdivide it into a finite number of small volume elements such 
as tetrahedra, triangular prisms, or rectangular bricks. By using the edge-based vector 
basis functions [25], the electric field E or magnetic field H is expanded as 

N v 

E (or H) = J2 E i (° r H i ) W i ( 12 ) 

j - 1 

where N v denotes the total number of element edges resulting from the subdivision 
including those on the surface S. Also, E 3 (or Hj ) denote the unknown expansion 
coefficients equal to element edge fields and W j are the chosen vector basis functions. 
Note that the same basis functions are employed for both E and H, though this is not 
required. 

To generate a system of equations for the fields in V, (12) is substituted into (5). 
Applying the Rayleigh-Ritz procedure to enforce (6) then yields 


[>1]n v xN v {-E}jv v xi + [B]n v xn 9 {Hs}n 3 xi = {C}/v v xi (13) 

where {F} = [Fi, Ejv v ] t , {Hs} = [Hu Hi , ...» Bn b ] T with the superscript T denot- 

ing the transpose of the vector and N s being the total number of element edges residing 
on the surface S. The elements of [ A ], [B] and { C } are given by 


*ij 


Bij 


Ci 


JJIv (V X W,) # (V X Wj) ~ k2 ° €rWi • Wj ] 

JkoZ 0 jj) s W.-CWj x n)dS 



jk 0 Z 0 J int - V x 



dV 


(14) 

(15) 

(16) 
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We note that the result stated in (13) can also be obtained via a discretization of the 
weighted residual equation (9) provided the same expansions are used and T is set to 
T = W„ implying a Galerkin’s formulation. Finally, a discretization of (10) or (11) 
yields the dual of (13) given by 

}jv„xi + [•B , ];v ) , X Ar 1 {£s}/'/ a xi = {C'JiVvxi • (17) 

The elements of the matrices [ A '], \B'] and {C'} are, respectively, given by the dual of 
(14)-(16). 

4. Boundary Element Discretization 

To solve the system (13) or (17), a specification of a boundary condition relating 
(E X n) and (H X n) is required. This is provided by the boundary integral equation 
(1), and to incorporate it into (13) or (17) it must first be discretized. To illustrate this 
let us rewrite (1) for the electric field as 

E = E*" c + L a el (E x n) + L* 2 (H x n) (18) 

where L® was split into two parts; one (L®j ) pertaining to the equivalent magnetic current 
and the other (L* 2 ) to the equivalent electric current. Taking the cross product of (18) 
with n yields 

E x n = E mc x n + Lj X (E x n) x n + L* 2 (H x n) x n (19) 

which can be discretized by using (12) and applying Galerkin’s procedure to find 

[^JjVjxTV, {Es}n,xi + [■f > ]N a x^ s {-£'s}^ J xi + [Q]n,xnAHs}n,xi = {^}tv iX i • (20) 


8 



The matrix elements of [ B * ] are the same as those in (17) and the other matrix elements 
are given by 


Pij = jk 0 Y 0 

W, • [L* a (Wj x n) x n] dS 
s 

(21) 

Qij = iWjg 

W i • [L *2 C Wj x n) x n] dS 

> 

(22) 

Y = -jk 0 Yo j 

fj)^ W, • (E ,nc x n) dS . 

(23) 


The dual of (20) can also be obtained from (1) and is given by 

[B]n.xN,{Hs}n,x 1 + \P']n.xN,{Hs}n s xi + [<5V»xiV,{^5}^xl = • (24) 

In this, the elements of [5] are the same as those in (13) and the elements of [P'], [Q'\ 
and [ Y '] are given by (21)-(23), respectively, upon replacing Yo with — Zq, U eX 2 with 
L£ 1>2 and E inc with H ,nc . 

5. Solution of the System 

A complete system of equations for the discrete edge fields can now be obtained by 
combining any one of (13) and (17) with either (20) or (24). Various techniques can be 
employed for solving the resulting system, but it would be advantageous to use algorithms 
which exploit the special properties of the finite element matrices. More importantly, the 
subsystems (20) and (24) resulting from the discretization of the boundary integrals are 
fully populated and could increase the memory demand beyond O(N) unless special care 
is exercised. In particular, to retain the O(N) storage requirement, the boundary/surface 
S must be judiciously chosen so that the resulting boundary integral or a large portion 
of it is convolutional. The FFT can then be used to evaluate the integral without a need 
to explicitly generate and store the boundary element matrices, provided an iterative 
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solution such as the CG or BiCG method is employed. Below we consider specific two- 
and three-dimensional applications which illustrate how the boundary /surface S could 
be chosen to maintain the 0(N ) storage requirement. These applications will also permit 
the derivation of specific forms for the matrices introduced in (13), (17), (20) and (24) as 
well as demonstrate the validity and effectiveness of the method. In the following we first 
consider a few two-dimensional applications and in this case V and represent cross- 
sectional areas denoted by Q. and Q , respectively, and S becomes a contour which will be 
denoted by I\ Some three-dimensional applications are also discussed and computations 
are presented which are compared with reference data. 

III. TWO-DIMENSIONAL ANALYSIS 

For two-dimensional problems, it is sufficient to consider the transverse magnetic (TM 
or i?-polarization) and transverse electric (TE or 77-polarization) incidence separately. 
This reduces the problem to a scalar one and (5) or (10) becomes 

m= \!L{ u (j £) + (§£) _ - * 2 *^} dxiy + 1 ** ir < 25 > 

where 

4> - E z , xp = jk 0 Zo(K x n) • z , « = — , v = c T (26) 

f*T 

for ^-polarization and 

4>= H z , xp = -jk 0 Yo( E x n) • z , u = — , v = n T (27) 

for /7-polarization. The matrix equations (13) and (17) now reduce to 

[A]{<P} + [B]{xP} = 0 (28) 
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where 




B , 




■ //, [• 
■ i 


dNf_dN[ dNf dN[ 
dx dx dy dy 


- klvNfN* 


dxdy 


l 9 x l] dr . 


(29) 

(30) 


where Nf and Xf denote the chosen basis functions for the interior and boundary fields, 
respectively, and furthermore, Nf becomes the same as Xf at the boundary T. It remains 
to obtain the elements associated with (20) or (24). Obviously, the explicit form of these 
depends on the formulation of the boundary integral equation which is dependent upon 
the geometry of the scatterer under consideration. 


1. Scattering by Grooves and Slots in a Thick Conducting Plane 

Consider the geometry illustrated in Figure 2 where an infinitely long groove is cut 
in an infinite ground plane. For this specific configuration, ft is the region occupied by 
the cross-section of the groove and ft ^ is the half space above the ground plane (y > 0). 
Also, T consists of the straight line segment T i across the opening and the conducting 
boundaries forming the groove. Because of the boundary condition, the portion of the 
boundary integral in (25) over the conducting boundaries vanishes and, thus, we only 
need to consider the remaining boundary integral along T \ . Below we discuss the TM and 
TE cases separately since they are associated with different boundary integral equations. 

A. TM incidence 

Using the notation defined in (26), the boundary integral equation is given by 

v>(x, 0) = 2ip inc (x,0) + J - ^*2 + 4>(x')H^ 2) {ko\x - x'\)dx' (31) 
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where denotes the zeroth order Hankel function of the second kind and a factor of 

two has been introduced in the right-hand side to account for the presence of the ground 
plane. To obtain a system corresponding to (24), (31) is discretized via Galerkin's 
technique, yielding 

= (32) 

In this, {fo} denotes the portion of {0} on the boundary T i, [5] is the same as that in 
(28) and [Q f ] and {Y*} are given by 

Q % = - 5 //?(*) + \x-x'\)dx' dx (33) 

Y( = 2 [ L 9 t (x)rp ,nc {x,0)dx (34) 

which can be evaluated analytically or numerically for the given Zf and L 9 . 

To solve for the fields {<£}, (32) and (28) can be combined to yield 

- mm = -<n (35) 

which is amenable to a unique solution and can be solved via a number of methods. 
However, if the CG or BiCG method is used, the FFT can be employed for the compu- 
tation of [Q f ]{4>b] without a need to generate the matrix [Q f ]. To illustrate this we refer 
to (33) and observe that Q is a function of (i — j) for an equal subdivision of T j . Thus 
we may invoke the convolution theorem to write 

IQ') m = ?-'{? W}oT{4> b }} (36) 

in which q\ = Q' ix , T denotes the discrete Fourier transform and the symbol o implies 
the Hadamard product. Clearly, the use of the FFT in (36) eliminates a need to store 
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the entire [Q f ] matrix (other than one row of it), thus maintaining the O(N) memory 
requirement. We now turn to the computation of [A]{<£}. As noted earlier, [A] is 
an N v X N v square matrix resulting from the finite element discretization. Since it 
is assembled from the element matrices, its associated matrix product can be easily 
computed as 

M v 

MW = £M]W> (37) 

e=l 

where M v is the total number of subdivision elements, [ A e ] is the element matrix and 
{0 C } denotes the fields associated with the eth element. Note that [A e ] are generally 
simple small matrices and thus the computation of (37) requires only 0(N) memory. 
The total memory demand for the solution of the system (35) is then kept to 0(N ) since 
[A]{<£} and [£?']{<&>} are the only computations in a CG or BiCG algorithm involving 
the use of [A] and [Q*], 

B. TE incidence 

The procedure for this excitation is very similar to that outlined for the TM case. 
The boundary integral equation now is 

4>{x,0) = 2<t> mc (x,0)+ ^ - x'\)dx 9 (38) 

where the notation defined in (27) has again been adopted. A discretization of this via 
Galerkin’s procedure then yields 

(39) 

where [5] is the same as that in (28) and the elements of [Q l ] and [Y f ] are now given by 
Qij = !?(*) [ / L]{x') H ^\k 0 \x - ^1)^'] dx (40) 
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Y/ = 2 f Lf(x)4> inc (x,0)dx . (41) 

Jr i 

A solution for {^} and then follows by combining (39) with (28). As in the TM 
case we again observe that [Q f ]{ip} can be computed via the FFT as 

[<3lW = r ] {;{ ? ')o;w) (42) 

where q f = Q' {1 . Provided the CG or BiCG method is used, the memory requirement of 
the solution is only O(N). 

C. Numerical results 

The above formulations for the TM and TE scattering by a groove and slot have been 
implemented using linear expansion functions. To show the capability of the method, 
we consider the plane wave scattering by two different structures. Figure 3 shows the 
backscatter radar cross section (RCS) for a rectangular groove, compared with measured 
data [26]. Figures 4 and 5 show the transmission coefficient as a function of frequency 
and bistatic scattering patterns for a structure consisting of a wide slot having a non- 
uniform filling and containing a periodic array of strips in a multilayer dielectric. We 
note that the model employed for generating the data in Figure 4 required nearly 2500 
unknowns at the high end of the spectrum. However, by using the CG or BiCG method 
in conjunction with the FFT, we were able to carry out the solution of the system on a 
workstation. 

2. Scattering by Cylinders 

We consider now a cylindrical scatterer of arbitrary cross-section enclosed by the 
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fictitious boundary T . The boundary integral equation then becomes 

4>{p) = 4> inc {p) + j( [oiP') 990 ^/ 1 + <7o(p,pW)] dV (43) 

where p and p' denote the observation and integration points, respectively, and 

<7o(p,p') = - J -H ( Q ] {k 0 \p- p'\) (44) 

is the free space Green’s function. A discretization of (43) via Galerkin’s approach yields 

mh) + [p]m + [ om = (n (45) 

where [5] is given by (30) and the matrix elements for [P], [Q] and [F] are given by 

7^^]^ (^) 

Qa = - / r mp) [/ r ^(p')po(p, p'Mr'] dr (47) 

Yi = f r L°(p)<t> inc (p)dT (48) 

with p denoting the observation point on T. 


The combined system of (45) and (28) now forms a complete system for the solution 
of {<f>} and {ip}. However, care must be exercised for the computation of [P]{<^ fc } and 
[Q]W to ensure an 0{N ) memory demand for the entire system solution if the CG 
or BiCG method is used. In the previous examples, this was achieved by exploiting 
the convolutionality of the integral operators. This property is strongly dependent on 
the choice for T. It was shown that for planar T the boundary integral operators are 
convolutional, and this also holds for circular boundaries as well. Choosing T to be a 
circle tightly enclosing the cylindrical scatterer and subdividing it equally, (46) and (47) 
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become 


Pij 


Qij 


A /* 2 tT r /‘ 27 T 

~4 k ° a2 l l ’ML W> 


■ H 


( 2koa [sin 

2ir 


<£JZ_!£_ 

2 

r r2ir 


sin 


<P-<P 


dtp' 


dtp 

<P-<? 


4°* Jo L *^'J ( 2 *° q sin ~ 2 ~ ) d(p ' 


dip 


(49) 

(50) 


where a denotes the radius of the circle and /Tj (•) denotes the first order Hankel 
function of the second kind. It can be readily shown that and Qij are functions 

of (i — j) and thus the matrix products associated with the boundary integrals can be 
efficiently evaluated via the FFT. Specifically, we have 


[P]{<t>b} = 

V) 

i 

r- A —\ 

Si 

l 

o 

Si 

& 

Cr* 

(51) 

[Q]M = 

o/W] 

(52) 


where pi = Pn and g, = Q tl . In passing, we should note, though, that choosing a circular 
boundary enclosure may not necessarily be a memory efficient approach, particularly if 
the scatterer is small in one dimension. In that case, a rectangular or ogival enclosure that 
tightly encloses the target will most likely result in less unknowns [20]. The boundary 
integrals are then convolutional only when the observation and integration points are on 
the same side of the ogive or the rectangle, and also when they are on the parallel sides 
of the rectangle. Otherwise, the boundary integral has no special form and its associated 
matrix must therefore be stored explicitly [16]. 

The above formulation has been implemented and validated for a variety of geome- 
tries. Figure 6 shows the backscatter patterns for a coated ogival cylinder in comparison 
with the moment method data [27]. In generating the finite element - boundary integral 
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solution we employed a ogival enclosure in order to keep the finite element region to 


minimum. 


IV. THREE-DIMENSIONAL ANALYSIS 

The full advantage of the proposed hybrid method is realized when we consider a 
three-dimensional application where the memory demands are far more excessive. There- 
fore, having a solution method leading to an 0(N ) memory demand is of crucial impor- 
tance. In this section, we consider the problem of scattering by a cavity-backed aperture 
and a slot in a thick conducting plane, scattering and radiation by a microstrip patch 
antenna or array in a cavity, and that of scattering by a finite size object. 

1. Scattering by Cavity-Backed Apertures and Slots in a Thick Conducting 
Plane 

Consider the cavity-backed aperture illustrated in Figure 7. In this case, is the 
free space region above the ground plane {z > 0) and V is that occupying the cavity 
(— c < z < 0). The surface S consists of the planar aperture and the conducting walls 
of the cavity. Because of the boundary condition the portion of the boundary integral 
over the conducting walls vanishes as was the case with the groove. Thus, we only need 
to consider the remaing portion of the integral over the aperture. 

For this problem, the boundary integral equation (1) is given by 

H(r) = H' nc (r) + H re '(r) - 2jk 0 Y 0 JJ G 0 (r,r') . [E(r') x z] dS' (53) 

where H* nc denotes the incident field due to (3 ext , M ext ) and H re/ is that reflected 
by the ground plane without the aperture. Also, G 0 is the free-space dyadic Green’s 
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function defined by 


G 0 (r,r')= 1 + 


G 0 (r,r') 


( 54 ) 


with 


Go(r,r') = 


e — J*0 |r— r'l 
47rir - r'l ' 


( 55 ) 


Multiplying (53) by (jk 0 Zo), crossing by z and discretizing the resulting integral equation 
via Galerkin’s method yields 


[£]{tfs} + [ Q']{Es } = {¥'} . (56) 

In this, [ B ] is given by (15), \Q f ] and {Y'} are similar to those in (24) and are more 
explicitly given by 

Q'ij = 2^o JJ s [Wi(r) X z] • \^jj s Co(r,r ; ) • [Wj(r') x z] dS ( 57 ) 
Y/ = 2jk 0 Z 0 JJ v/i(r)* [H ,nc (r)xz] dS . ( 58 ) 

Obviously, the integrand singularty of Q • is nonintegrable and it is necessary to employ 
the divergence theorem in order to transfer the del operators contained in Go to the 
expansion and weighting functions and this is discussed in [21]. As before, a final system 
of equations is obtained by combining (56) and the finite element equation (13) to give 

M{£}-[Q']{£s} = -{n- (59) 

This can be solved using various techniques including the CG or BiCG method in con- 
junction with the FFT for the evaluation of the product [< Q 'J {£s}- 
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The above formulation was implemented and validated using linear basis functions for 
rectangular cavities and slots which are amenable to simple finite element discretizations. 
As an example, Figure 8 displays the backscatter RCS of a rectangular cavity and as 
seen, the results based on this formulation are in excellent agreement with measured 
data [28]. Other computations are displayed in Figures 9 and 10. In particular, Figure 
9 presents the scattering by a circular cavity whereas Figure 10 displays the bistatic 
scattering by a material filled rectangular slot. Of course, the presented formulation is 
applicable to cavities filled with inhomogeneous material whereas traditional approaches 
are not and one such application is considered next. 

2. Scattering and Radiation by Microstrip Patch Antennas in a Cavity 

The structure to be considered is illustrated in Figure 11 where a microstrip patch 
antenna or array is residing on or embedded in a substrate which is in turn housed in 
a cavity recessed in a ground plane. As expected, the formulation for this problem is 
similar to that of scattering by a cavity-backed aperture, except in the case of radiation 
where the excitation is due to internal sources in the cavity. The system to be solved is 
therefore 

[A]{E}-[Q']{E S } = {C} (60) 

rather than (59) which is suitable for scattering computations. In (60), {C} is given 
by (16) and the product [Q'] {£s} is again evaluated via the FFT for a CG or BiCG 
solution. 

The modeling of the conducting patches and microstrip transmission lines is carried 
out by setting the electric field components to zero for those element edges coinciding with 
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the patch. Further, an impedance load can be modeled as a post of finite conductivity 
joining both the patch and the base of the cavity. The antenna can be excited either by 
a current filament, a magnetic frill current, or a gap generator [23]. 

Figure 12 shows the backscatter RCS spectrum for a single patch loaded with 50 
ohms at the point (xl>2/l)- About 1280 unknowns were required for the simulation 
of the geometry. Notably, on the average only 100 iterations were required for the 
solution to converge to within 0.01 dB of the correct RCS level for each excitation with 
a corresponding cpu time of about 20 s on a Cray Y-MP832. The radiation patterns for 
a 13x16 microstrip patch array is given in Figure 13. The patches of this array are of 
the same size as that shown in Figure 12. They are 1.83 cm apart in the x-direction, 
1.30 cm apart in the y-direction, and the cavity dimensions are 73.2 cm x 63.7 cm x 
0.158 cm. Each patch is uniformly fed at its lower left corner and the feed was modeled 
by a current filament. For this example, 120935 unknowns were used and the solution 
converged within 100 iterations. 

3. Scattering by a Finite Body 

The corresponding equation to (1) for a finite body enclosed by the fictitious surface 
5 is given by 

E(r) = E ,nc (r) + {v x G 0 (r, r') • [E(r') x n'] 

— jk 0 ZoOo(*, 1*0 • [H(r') X n']} dS' (61) 

for the electric field. A numerical discretization of (61) yields the matrix equation (20), 
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whose matrix elements can be obtained from (21)-(23). Specifically, we have 

Pij = -J k oYofy f [W,(r) X n] . V x G 0 (r, r') • [W.(r') x n'] dS'} dS (62) 
Qij = -*o j| [W,(r) x n] . {j| G 0 (r,r') . [W,(r') x n'] d5'} <i5 (63) 

The implementation of the above formulation for an arbitrary body has not yet been 
carried out. A major difficulty is the large number of unknowns on S even for small 
bodies. To overcome this, a cylindrical enclosure may be chosen so that the integrals over 
the cylindrical surface become two-dimensional convolutions and can, thus, be evaluated 
via the FFT, provided the CG or BiCG method is employed for the solution of the 
system. This idea has been tested for a body of revolution [24] and Figure 14 displays 
the bistatic scattering patterns for an ogive as obtained by this method and the moment 
method [29]. 

V. CONCLUSION 

In this paper we reviewed a hybrid finite element - boundary integral method for 
scattering and radiation applications. The method involves the introduction of a ficti- 
tious boundary enclosing the scatterer which serves to decouple the fields in the regions 
interior and exterior to the boundary. The fields in the interior region are formulated 
via the finite element method whereas those in the exterior region are represented by the 
boundary integrals involving the free space Green’s function. The interior and exterior 
fields are then coupled by invoking the continuity of the tangential boundary fields, re- 
sulting in a complete system for the solution of the fields internal to and on the fictitious 
boundary. 

Of particular interest in the proposed hybrid formulation is the choice of the fictitious 
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boundary so that the resulting boundary integrals are convolutional. If so, they can 
then be evaluated via the FFT when an iterative solution of the system such as the 
conjugate gradient or biconjugate gradient method is employed. This avoids a need 
to generate the matrix and since the finite element discretization results in sparse and 
banded matrices, the entire system solution requires alow 0(N ) storage. The proposed 
technique, therefore, holds a promise for treating large structures without a compromise 
in accuracy. We should note, however, that this technique, like others involving the use 
of boundary integral equations over closed surfaces or contours, may be associated with 
fictitious internal resonance phenomena which can be eliminated by combining the E 
and H field boundary integral equations. Such difficulties do not, of course, arise when 
the boundaries are not closed such as in the case of aperture problems. 

Finally, we note that in all cases considered in this paper the final system of equations 
is symmetric, due to the employment of Galerkin’s method for the discretization of the 
boundary integral equations. As a result, the BiCG method is more favorable than 
the CG method since it requires only one matrix-vector product computation in each 
iteration and also since it converges faster. 
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FIGURE CAPTIONS 


Fig. 1 Geometry of the scattering/radiation problem. 

Fig. 2 Geometry of a two-dimensional groove. 

Fig. 3 Backscatter RCS for a 2.5 cm wide and 1.25 cm deep rectangular groove as a 
function of frequency, if-polarization, y? tnc = 10°. 

Fig. 4 Transmission coefficient as a function of frequency for a truncated strip grating 
at normal incidence. Top and bottom layers: € r = 2.56, ji r = 1.0, 0.2 cm thick; 
middle layer: e r = 4.0, \x T = 1.0, 0.2 cm thick. 

Fig. 5 Bistatic scattering patterns for the geometry in Figure 4 at 10 GHz and y> xnc = 
60°. (a) F-polarization; (b) if-polarization. 

Fig. 6 Backscatter pattern for a 4A x 1A perfectly conducting ogival cylinder with a 
0.05A thick material coating having e r = 3 — j5,/x r — 1.5— j'0.5. (a) F-polarization; 
(b) if-polarization. 

Fig. 7 Geometry of a cavity-backed aperture in a ground plane. 

Fig. 8 Backscatter pattern for a 16 inch long, 0.1968 inch wide and 0.837 inch deep 
cavity at / = 12 GHz. E ,nc = y F, </? inc = 0. 

Fig. 9 Backscatter pattern for an 1 inch deep circular cavity having a diameter of 1 
inch at / = 16 GHz. 
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Fig. 10 Bistatic scattering pattern at normal incidence for a 1.2A long, 0.25A wide and 
0.25A deep slot filled with the material of e r = 2 — j and (i r = 1.2 — j’0.1. E tnc = x£ 
and circles correspond to the scaled two-dimensional solution [19]. 

Fig. 11 Geometry of a microstrip patch array in a cavity. 

Fig. 12 Backscatter RCS versus frequency for a single patch loaded with 50 ohms. The 
cavity is 0.158 cm deep and filled w'ith a substrate having e r = 2.17 and a loss 
tangent of 0.001. 9 inc = 60°, 9* nc = 45°, E ,ne = BE. 

Fig. 13 Radiation pattern of the 13 x 16 microstrip patch array equally fed at the lower 

left corner of each patch. ( — ) <p = 0 — 7r plane; ( ) (p = 7r/2 — 37t/2 plane, (a) 

/ = 2.62 GHz. (b) / = 3.55 GHz. 

Fig. 14 Bistatic scattering pattern for a 1.0A long conducting ogive having a diameter 
of 0.176A at its center. 
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